Magnetoinductive breathers in magnetic metamaterials 
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The existence and stability of discrete breathers (DBs) in one-dimensional and two-dimensional 
magnetic metamaterials (MMs), which consist of periodic arrangements (arrays) of split-ring res- 
onators (SRRs) , is investigated numerically. We consider different configurations of the SRR arrays, 
which are related to the relative orientation of the SRRs in the MM, both in one and two spa- 
tial dimensions. In the latter case we also consider anisotropic MMs. Using standard numerical 
methods we construct several types of linearly stable breather excitations both in Hamiltonian and 
dissipative MMs (dissipative breathers). The study of stability in both cases is performed using 
standard Floquet analysis. In both cases we found that the increase of dimensionality from one 
to two spatial dimensions does not destroy the DBs, which may also exist in the case of moderate 
anisotropy (in two dimensions). In dissipative MMs, the dynamics is governed by a power balance 
between the mainly Ohmic dissipation and driving by an alternating magnetic field. In that case it 
is demonstrated that DB excitation locally alters the magnetic response of MMs from paramagnetic 
to diamagnetic. Moreover, when the frequency of the applied field approaches the SRR resonance 
frequency, the magnetic response of the MM in the region of the DB excitation may even become 
negative (extreme diamagnetic). 

PACS numbers: 63.20.Pw, 75.30.Kz, 78.20.Ci 
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I. INTRODUCTION 



Discrete breathers (DBs), also known as intrinsic lo- 
calized modes (ILMs), are genuine nonlinear excitations 
that oscillate for long times in a localized region of space 
and may be produced generically in discrete lattices of 
weakly coupled nonlinear elements (see [1| for a general 
review). Since their introduction |2j, a large volume of 
analytical and numerical studies have explored the exis- 
tence and the properties of DBs in a variety of discrete 
nonlinear systems. Rigorous mathematical proofs of ex- 
istence of DBs in both Hamiltonian (i.e., energy con- 
served) and dissipative lattices of weakly coupled non- 
linear oscillators have been given 0, and numerical 
algorithms for their accurate construction have been de- 
signed @, 0, 0|- DBs may appear spontaneously in a 
lattice as result of fluctuations [1, Tiof. disorder [TTI ] , 
or by purely deterministic mechanisms They have 

been observed experimentally in several systems, includ- 
ing solid state mixed-valence transition metal complexes 
[13| , quasi-one dimensional antiferromagnetic chains 14 1 , 
arrays of Josephson junctions [15J, micromechanical os- 
cillators [T(|, optical waveguide systems 17 1, and pro- 
teins (l8| . Once generated, DBs modify system proper- 
ties such as lattice thermodynamics and introduce the 
possibility of nondispersive energy transport [H, H3] > be- 
cause of their potential for translatory motion (i.e., mo- 
bility) along the lattice [2lj . In numerical experiments 
DB mobility can be achieved by an appropriate pertur- 
bation f22| ]. Recently, it has been found experimental 



evidence for moving DBs in a layered crystal insulator 
at 300-ftT [23| and in a small micromechanical cantilever 
array 24]. From the perspective of applications to ex- 
perimental situations where dissipative effects are always 
present, dissipative DB excitations (usually driven by an 
alternating power source) are more relevant than their 
Hamiltonian counterparts. Dissipative DBs, which pos- 
sess the character of an attractor for initial conditions in 
the corresponding basin of attraction, may appear when- 
ever power balance, instead of energy conservation, gov- 
erns the nonlinear dynamics of the lattice. Furthermore, 
the attractor character of dissipative DBs allows for the 
existence of quasiperiodic and even chaotic DBs [2|,[2^. 

Recently, dissipative DBs have been demonstrated nu- 
merically in discrete and nonlinear magnetic metamateri- 
als (MMs) driven by an alternating electromagnetic (EM") 
field [27j • The MMs are artificially structured materials 
that exhibit EM properties not available in naturally oc- 
curring materials. The response of any material (cither 
natural or artificial) to applied EM fields is characterized 
by macroscopic parameters such as the electric permit- 
tivity e and the magnetic permeability /x. For example, 
there are only a few natural materials responding mag- 
netically at Terahertz (THz) and optical frequencies, and 
that response is usually very weak. However the MMs 
exhibit relatively large magnetic response at those fre- 
quencies [111 [29|, [3fj| . which may be either positive or 
negative, resulting in positive or negative /x, respectively. 
The realization of MMs at such frequencies will certainly 
affect THz optics and its applications, while it promises 
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new device applications. The most common element uti- 
lized for magnetic response from MMs is the split-ring 
resonator (SRR) which, in its simplest form, is a metallic 
ring with one slit, made of a highly conducting metal. 
Periodic SRR arrays in the nanoscale, a genuine realiza- 
tion of MMs, are fabricated routinely using conventional 
microfabrication techniques. These MMs, in which the 
SRRs are weakly coupled magnetically through their mu- 
tual inductances, support a new t ype of guid ed waves, 
the magnetoinductive (MI) waves [3l|, HIH [Hj]. The 
MI waves propagate within a band near the resonant fre- 
quency of the SRRs, and exist as forward and backward 
waves depending on the orientation of the elements (e.g., 
the SRRs) of the MM. In the linear regime of MI wave 
propagation in a MM, the magnetic permeability fj, does 
not depend on the intensity of the EM field. However, 
the MMs may become nonlinear, either by embedding 
the SRRs in a Kerr-type medium [3^, [36[ , or by insert- 
ing certain nonlinear elements (e.g., diodes) in each SRR 
[37I l38l l39l ]. Then, the combined effects of nonlinear- 
ity and discreteness (inherent in SRR-based MMs) , leads 
in the generation of nonlinear excitations in the form of 
DBs (23] , as well as magnetic domain walls [4(| , and mag- 
netoinductive envelope solitons [4l|. The latter may be 
both bright or dark, and they result from the modula- 
tional instability of the MI wares. While the nonlincar- 
ity by itself offers tunability, stationary (i.e., not mobile, 
pinned) DBs act as stable impurity modes that are dy- 
namically generated and may alter propagation and emis- 
sion properties of a system. Moreover, stationary dissi- 
pative DBs can change locally the magnetic response of a 
nonlinear MM (27[. We should also note that regular ar- 
rays of rf SQUIDs offer an alternative for the construction 
of nonlinear MMs due to the nonlinearity of the Joseph- 
son junction (42T |. 

In the present work we investigate numerically the ex- 
istence and stability of both dissipative and energy con- 
served DBs in discrete, periodic arrays of nonlinear SRRs. 
We consider different SRR array geometries (i.e., differ- 
ent orientations of the SRRs in the MM) in one and two 
spatial dimensions. In two dimensions we also consider 
several cases of possible anisotropy. In the next section 
we describe the discrete MM model, while in section III 
we construct several types of DBs for one-dimensional 
(ID) arrays. Here we also calculate the magnetic re- 
sponse, which may be locally altered by the presence of a 
DB. In section IV we construct several types of DBs for 
two-dimensional (2D) arrays, and we finish in section V 
with the conclusions. 



II. MAGNETIC METAMATERIAL MODEL 

We consider ID and 2D discrete, periodic arrays of 
identical nonlinear SRRs, which consist the simplest re- 
alization of a MM in one and two dimensions [43[ , respec- 
tively. In one dimension, the SRRs form a linear array 
with their centers separated by distance d. In two dimen- 



sions, the SRRs are assumed to be arranged in a regular 
rectangular lattice, with their centers separated by dis- 
tance d x and d y in the x— and y— directions, respectively 
(i.e., lattice constants d x and d v , respectively). There 
are two configurations of interest in each case, shown 
schematically in Figs. 1 and 2 for the ID and 2D ar- 
rays, respectively. A ID array can be constructed either 
in the planar configuration, where all SRR loops are in 
the same plane with their centers lying on a straight line 
(Fig. 1(a)), or in the axial configuration, where the line 
connecting the centers of the SRR loops is perpendicu- 
lar to the plane of the loops (Fig. 1(b)). Similarly, a 
2D array can be constructed either in the planar config- 
uration, where all SRR loops are in the same plane with 
their centers located on an orthogonal lattice (Fig. 2(a)), 
or in the planar-axial configuration, where all SRRs have 
the planar configuration in one direction while they have 
the axial configuration in the other direction (Fig. 2(b)). 

Ay (a) 
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FIG. 1: Schematic view of a one-dimensional array of split- 
ring resonators in (a) the planar geometry; (b) the axial geom- 
etry. In both geometries the split-ring resonator axes as well 
as the magnetic component of the applied field are directed 
along the y— axis. The electric field component is transversal 
to the slits (not shown in the figure). 

Within good approximation, each SRR is equivalent 
to a nonlinear resistor-inductor-capacitor (RLC) circuit 
featuring a self-inductance L, Ohmic resistance R, and 
capacitance C . The units (i.e., the SRRs) become nonlin- 
ear due to a Kerr dielectric filling the SRRs' slits, whose 
permittivity e is of the form 

£ (|E| 2 )=e ^ + a^, (1) 

where E is the electric field, E c is a characteristic (large) 
electric field, q is the linear permittivity, eo is the per- 
mittivity of the vacuum, and a = +1 (—1) correspond- 
ing to self-focusing (self-defocusing) nonlinearity. As a 
result, the SRRs acquire a field-dependent capacitance 
C(|E| 2 ) = e(|E g | 2 ) A/dg, where A is the area of the cross- 
section of the SRR wire, E g is the electric field induced 
along the SRR slit, and d g is the size of the slit. The ori- 
gin of E g may be due to the magnetic and/or the electric 
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components of the applied EM field, depending on the 
relative orientation of that field with respect to the SRRs' 
plane and the slits. In the following we assume that, for 
any array configuration and number of dimensions, the 
magnetic component of the incident (applied) EM field 
is always perpendicular to the SRRs' plane, and that the 
electric component of the incident EM field is transversal 
to the slit. Then, only the magnetic component excites 
an electromotive force (emf) in the SRRs, resulting in 
an oscillating current in each SRR loop and the corre- 
sponding development of an oscillating voltage difference 
U across the slits or, equivalently, of an oscillating elec- 
tric field E ff in the slits. If Q is the charge stored in the 
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FIG. 2: Schematic view of a two-dimensional array of split- 
ring resonators in (a) the planar geometry; (b) the planar- 
axial geometry. In both geometries the magnetic component 
of the applied field is directed along the SRR axes, while the 
electric field component is transversal to their slits (not shown 
in the figure). 

capacitor of an SRR then, from the general relation of a 
voltage-dependent capacitance, C(U) — dQ/dU, and Eq. 
(P), we get 



C, 1 



u 2 



U, 



(2) 



where U — d g E g , Ce = eoei(A/d g ) is the linear capaci- 
tance, and U c = d g E c . Assume that the arrays are placed 
in a time- varying and spatially uniform magnetic field of 
the form 



H = H o cos(wf), 



(3) 



where Hq is the field amplitude, u> is the field frequency, 
and t is the time variable. The excited emf £ , which is 
the same in all SRRs, is given by 



£ = So sin(wi), £ a — A*o w S H( h 



(4) 



neighbor interactions, so that neighboring SRRs are cou- 
pled through their mutual inductances M x and M y . This 
is a very good hypothesis in the planar configurations 
(i.e., both in ID and 2D arrays), even if the SRRs are 
very close. The validity of the nearest neighbor approxi- 
mation for the other configurations (i.e., the axial config- 
uration in ID arrays and the planar-axial configuration in 
2D arrays) has been checked by taking into account the 
interaction of the SRRs with their four nearest neighbors, 

(s) 

using that the mutual inductance M x . y between an SRR 

and its s— th neighbor goes like Mx\y — M x ^ y /s 3 . We 
found that for weak coupling, as we consider here, the 
results are practically the same with those obtained with 
the nearest neighbor approximation. Thus, the electri- 
cal equivalent of an SRR array in an alternating mag- 
netic field is that of nonlinear RLC oscillators coupled 
with their nearest neighbors through their mutual induc- 
tances, which arc driven by identical alternating voltage 
sources. Therefore the equations describing the dynam- 
ics of the charge Q n ,m and the current I n ,m circulating 
in the n, m— th SRR may be derived simply from Kirch- 
hoff's voltage law for each SRR [13, lioj 
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where f(Q n ,m) — U n _ m is given implicitly from Eq. 
Using the relations 

ujg 2 = Ld, t = twi, I c = U c Lu e Ce, Q c = CtU c (7) 

£ — UqS : In,m — ^c^n,m, Qn,m — QcQn.nii (8) 

and Eq. J!]), Eqs. (0) and (J6j) can be normalized to 

d>Qn rn 



dr 

din,m . »/ \ i \ / l,m ^n-f-l,m 

+ 7*71.771 + J\Qn,m) + *x 
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£ sin(Or), (10) 



where 5* is the area of each SRR loop, and /^o the per- 
mittivity of the vacuum. Each SRR in the field given by 
Eq. ([3]) is a nonlinear oscillator which exhibits a reso- 
nant magnetic response (either positive or negative) at 
a particular frequency which is very close to its linear 
resonance frequency ui£ — 1 / \JL Ci (for R ~ 0). 

All SRRs in an array are coupled together due to 
magnetic dipole-dipole interaction through their mutual 
inductances. However, we assume below only nearest 



where 
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RCiuji is the loss coefficient, X x 



M x>y /L are the the coupling parameters in the x— and 
y— direction, respectively, and £o = £q/U c . Note that the 
loss coefficient 7, which is usually very small (7 <C 1) , 
may account both for Ohmic and radiative losses [4JJ. 
The corresponding equations for ID arrays result from 
Eqs. (HJ) and (jTUJ) , i.e., by setting X y = 0, dropping the 
subscript m, and choosing the appropriate X x = X. Ne- 
glecting losses and without applied field, Eqs. ([5]) and 
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(IIOP can be derived from the Hamiltonian 




^ ^ {^x Qn.jn (?n+l,m 4" Ay Qn,m Qn,m+l} ; (-H) 

where the nonlinear on-site potential V n ,m is given by 

V n>m = V(q n , m ) = / f(q'n,m) d q'n,m, ( 12 ) 

Jo 

and <7„. m = dq nm /dr. Analytical solution of Eq. @ for 
w n ,m = f(Qn,m) with the conditions of u„ im being real 
and u„ im (<?«., m = 0) = 0, gives the approximate expres- 
sion 

f(q n , m ) ^ 9n,m ~ 3^«n,m + 3 f 3^ J ?n,m> ( 13 ) 

which is valid for relatively low q n (q n < 1, 71 = 
1,2,..., AT) Thus, the on-site potential is soft for a — +1 
and hard for a = — 1. In the 2D case the mutual in- 
ductances M x and M y may differ both in their sign, de- 
pending on the configuration, and their magnitude. Ac- 
tually, even in the planar 2D configuration with d x = d y 
a small anisotropy should be expected because we con- 
sidered SRRs having only one slit. This anisotropy can 
be effectively taken into account by considering slightly 
different coupling parameters \ x and \ y . The coupling 
parameters \ x _ y as well as the loss coefficient 7 can be 
calculated numerically for this specific model within high 
accuracy. However, for our purposes, it is sufficient to es- 
timate these parameters for realistic (experimental) ar- 
ray parameters, ignoring the nonlinearity of the SRRs 
and the effects due to the (weak) coupling. The self- 
inductance of a circular SRR of radius a (with circular 
cross-section of diameter h) can be determined by the 
well-known expression L — /xoa[ln(16a//i) — 1.75]. The 
value of the capacitance then follows from the choice of 
the resonance frequency(~ Lug). The resistance can be 
calculated from the given SRR dimensions and the ap- 
propriate value of the conductivity, taking into account 
the skin effect. The expression for the mutual inductance 
between two loops can be calculated by means of a simple 
approximation [271 ] . For example, for an array of circular 
silver-made SRRs with circular cross-section in the pla- 
nar geometry, with geometrical and material parameters 
close to those in Ref. [43| (for the equivalent squared 
SRR with squared cross-section) . it has been estimated 
that A w 0.02 and 7 « 0.01 0. For the same SRRs 
in the axial geometry, separated by distance d as those 
in the planar geometry, the approximate expressions for 
the mutual inductance and the coupling parameter are 
M ~ (Tr/2)(j, a(a/d) 3 and A ~ (n/2)^ a(a/d) 3 /L, re- 
spectively. Those expressions are obtained by making 
the same approximations as in Ref. [27], and knowing 
that the magnetic field of one of the SRRs with induced 



current / at the center of the other SRR, which is directed 
along its axis, is given by B = ^qIo 2 /2(a 2 + d 2 ) 3 / 2 . Thus, 
in a first approximation, the coupling of two SRRs in the 
axial geometry is stronger (by a factor of 2) than that of 
the same SRF° irl fV, ° n1 '""" ™ nmotT " 
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FIG. 3: (Color online) The frequency spectrum of linear MI 
waves Q K as a function of the wavenumber k, for A = —0.01 
(black-solid curve) and A = +0.01 (red-dashed curve). 

For the integration of Eqs. (|9|) and (fT0|) . both in the 
Hamiltonian and the damped driven cases, or the cor- 
responding linearized ones, we use a standard fourth- 
order Runge-Kutta algorithm with fixed time-stepping 
At (typically At — 0.01). Since the DBs studied here are 
highly localized, the choice of boundary conditions to be 
imposed on Eqs. and (flTj|) is not especially impor- 
tant. Thus, we have chosen periodic boundary conditions 
throughout the study. 

III. DISCRETE BREATHERS IN 
ONE-DIMENSIONAL SRR ARRAYS 

We consider the two different ID configurations of SRR 
arrays shown in Fig. 1, with the same number of SRR 
oscillators N — 50. Within the equivalent circuit model, 
the difference between the two configurations resides in 
the sign and the magnitude of the coupling parameter 
A. In the planar geometry (configuration) A is negative, 
since the mutual inductance M between two neighbor- 
ing SRRs is negative. This is due to the fact that the 
magnetic field generated by one SRR (due to the current 
induced in its loop) crosses the neighboring SRR in the 
opposite direction. For the axial geometry, the mutual 
inductance M, and hence A is positive. Moreover, in a 
ID array of SRRs in the axial geometry, the value of A 
is much higher than that in a ID array of SRRs in the 
planar geometry with the same SRR spacing d. High 
coupling between SRRs in the axial geometry would pos- 
sibly require to take into account the interaction of an 
SRR with its far neighbors (and not only with its near- 
est neighbors). To avoid such complications, and for the 
sake of comparison of the DBs obtained in the two ge- 
ometries, we assume that the SRR spacing is such that 
the magnitude of the coupling between neighboring SRRs 
is the same (or at least of the same order) in both geome- 
tries. 

For Hamiltonian systems DBs may be constructed 
from the anti-continuous limit [5j , where all the SRRs are 
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uncoupled (A — > 0), obeying identical dynamical equa- 
tions. Fixing the amplitude of one of them (say the one 
located at n = rib) to a specific value qb, with the corre- 
sponding current % — 0, we can determine, either ana- 
lytically (if possible) or numerically the period of oscilla- 
tion Tb- An initial condition with q n = for any n ^ rib, 
Qn b = lb, and q n — i n — for any n, represents a triv- 
ial DB with period Tb- Continuation of this solution for 
A 7^ using the Newton's method @, results in numeri- 
cally exact DBs up to a maximum value of the coupling 
parameter A = X max . For the existence of Hamiltonian 
DBs it is required that the DB frequency ujb — 2n/Tb, as 
well as all its multiples, lie outside the linear dispersion 
band of MI waves. The MI wave band is obtained by sub- 
stituting a plane wave of the form q n — Acos(ku — fir) 
into the linearized equations Eqs. © and (fit))) (with 
e = 0,7 = 0) 



n« = 



i 



yJl + 2X cos(k) ' 



(14) 



where fl = to /u>£ is the normalized frequency and k — kd 
is the normalized wavenumber (— n < k < n). Typical 
dispersion curves are shown in Fig. 3 for both geome- 
tries. For |A| -C 1 the band is very narrow (bandwidth 
Aft ~ 2|A|), so that the requirement for the existence 
of Hamiltonian DBs can be easily satisfied. Clearly, the 
bandwidth increases with increasing A. The MI waves 
are forward in the axial configuration (A > 0) with co- 
directional phase and group velocities, and backward in 
the planar configuration (A < 0), with phase and group 
velocities in the opposite directions. 

Following the procedure described above, with the ap- 
propriate choice of initial conditions (trivial DBs), we 
have constructed several types of Hamiltonian, numeri- 
cally exact DBs for ID SRR arrays in both the planar 
and the axial geometries, shown in Figs. 4 and 5 for 
self-focusing (a = +1) and self-defocusing nonlinearities 
(a = —1), respectively. Those Hamiltonian DB profiles 
(shown at maximum amplitude), i.e., the normalized cur- 
rent i n as a function of array site n, are characterized as 
single-site bright DBs (Figs. 4(a) and 5(a)), antisymmet- 
ric bright DBs (Figs. 4(b) and 5(b)), single-site dark DBs 
(Figs. 4(c) and 5(c)), and bright multibreather (Figs. 
4(d) and 5(d)). The term "bright" ("dark") DBs is used 
when there are only one or a few SRRs in which the 
current oscillates with large (small) amplitude, whereas 
the rest of them oscillate with small (large) amplitude 
[44| . All these DBs are highly localized, occupying only 
a few lattice sites, since they have been obtained for 
low values of the coupling parameter A. They are all 
symmetric, except those in Figs. 4(b) and 5(b) which 
are anti-symmetric. It is interesting to observe that, for 
self- focusing nonlincarity (a = +1), the profile of the 
single-site bright DB (Fig. 4(a)) is staggered (unstag- 
gered) for SRR arrays in the planar (axial) geometry, 
while, for self-defocusing nonlinearity (a — —1), the pro- 
file of the single-site bright DB (Fig. 5(a)) is staggered 
(unstaggered) for SRR arrays in the axial (planar) geom- 
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FIG. 4: (Color online) Several Hamiltonian discrete breather 
profiles (at maximum amplitude) for a — + 1, fit = 0.938, 
ei — 2, N = 50, constructed with Newton's method for both 
the planar and the axial array configurations (shown as black 
circles and red squares, respectively.) (a) A single-site bright 
breather for |A| = 0.02; (b) an antisymmetric bright breather 
for |A| = 0.02; (c) a single-site dark breather for |A| = 0.002; 
and (d) a bright multibreather for |A| = 0.013. Only part of 
the simulated array is shown for clarity. 



etry. Recall that a staggered (unstaggered) DB is that 
whose profile is of the form i n — X n exp[inb(n — Ji&)], 
In = \in\ > 0, with Kb = 7T (Kb = 0) [451 ] - The single-site 
DBs in Fig. 4(a) and 5(a), with (normalized) frequencies 
fib = L^b/i^e = 0.938 and 1.056, respectively, may be con- 
tinued for higher values of coupling. They cease to exist 
when the MI wave band, which expands with increas- 
ing A reaches the DB frequency fib- That will occur at 
|A| = \X max \ = |1 - l/ngl/2, which gives \X max \ ~ 0.068 
and 0.052 for Q b = 0.938 and 1.056, respectively. 

The linear stability of Hamiltonian DBs is addressed 
through the eigenvalues of the Floquet matrix (Floquet 
multipliers). A DB is linearly stable when all its Floquet 
multipliers lie on a circle of radius unity in the complex 
plane. The DBs shown in Figs. 4(a), 4(b) and 5(a), 5(b) 
are all linearly stable. However, the dark DBs shown in 
Figs. 4(c) and 5(c) as red squares, as well as the multi- 
breathers shown in Figs. 4(d) and 5(d) as red squares and 
black circles, respectively, are linearly unstable. Those 
DBs were found to be linearly stable only for very low 
couplings. 

Next, we construct DBs for SRR arrays which are sub- 
jected to losses (damping) and the external driving force 
e(r) (dissipative DBs). In order to generate DBs in this 
case we start by solving Eqs. and (TTTJ|) in the anti- 
continuous limit, where all SRRs are uncoupled. We 
identify two coexisting (and stable) attractors of a single 
damped-driven SRR oscillator with focusing nonlincar- 
ity (a = +1) and Q = 0.92, e = 0.04, 7 = 0.01, which 
have different amplitudes q^ ~ 1.6086 and qi ~ 0.2866 
(high and low amplitude attractors, respectively). No- 
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FIG. 5: (Color online) Several Hamiltonian discrete breather 
profiles (at maximum amplitude) for a — —I, fij, = 1.056, 
eg — 2, N — 50, constructed with Newton's method for both 
the planar and the axial array configurations (shown as black 
circles and red squares, respectively.) (a) A single-site bright 
breather for |A| = 0.02; (b) an antisymmetric bright breather 
for |A| = 0.02; (c) a single-site dark breather for |A| = 0.001; 
and (d) a bright multibreather for |A| = 0.013. Only part of 
the simulated array is shown for clarity. 
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tice that the frequency f2 is below the SRR resonance 
frequency, where a MM is expected to show only posi- 
tive magnetic response (i.e., response corresponding to 
positive /i). Subsequently, we fix the amplitude of one 
of the SRR oscillators (say the one at n = n^) to qh 
and all the others to qe (i n are all set to zero) . Using this 
configuration (trivial dissipative DB) as initial condition, 
we start integrating Eqs. and (TlU|) . while increasing 
the coupling parameter A in small steps. In this way, 
the initial condition can be continued for A ^ lead- 
ing to dissipative DB formation The time evolution 
of typical dissipative (single-site) bright DBs is shown 
for SRR arrays in the planar and the axial geometries 
in Figs. 6 (top and bottom panels, respectively). We 
see that both the DB and the background are oscillating 
with different amplitudes (high and low, respectively). In 
this aspect, dissipative DBs differ from their Hamiltonian 
single-site counterparts, where the background is always 
zero. We should note here that the DB frequency in this 
case has to be identical with that of the driving field, so 
that fih = 2ir/Ti, = f2. However, the phase differences of 
the SRR oscillators in an array with respect to the driv- 
ing field are generally different, so that a DB may change 
locally the magnetic response of an array, as we shall see 
below. With appropriate initial conditions we can also 
obtain multi-breathers where two or more sites oscillate 
with high (low) amplitude, while the other ones with low 
(high) amplitude. 

The linear stability of the dissipative DBs can be ad- 
dressed (as in the Hamiltonian case) through the eigen- 
values of the Floquet matrix (Floquet multipliers). In 



FIG. 6: (Color online) Time evolution of a single-site bright 
dissipative breather during approximately two periods, for 
fl b = 0.92, £ = 0.04, 7 = 0.01, u = 2, N = 50, and (top 
panel) planar SRR array geometry with A = —0.02; (bottom 
panel) axial SRR array geometry with A = 0.017. 



the dissipative case, however, a DB is linearly stable 
when all its Floquet multipliers lie on a circle of radius 
R = cxp(—jTb/2) in the complex plane p, due to the 
presence of a dissipative term in the linearized equations 
of motion. Both the dissipative DBs shown in Fig. 6 are 
linearly stable. In order to check that result, we have 
also added small perturbations to these DBs (of the or- 
der of 10 -2 of the DB amplitude) and let them evolve in 
time. We have followed the perturbed DBs over 10 3 Tf, 
time units observing that the DBs restore their unper- 
turbed shape. In general, dissipative DBs have been 
found to exist for couplings A much less than those of 
their Hamiltonian counterparts, X m ax- The DBs shown 
in Fig. 6 for dissipative and forced SRR arrays in the 
planar and the axial geometries, are found to exist up to 
A ~ 0.024 and 0.017, respectively. For defocusing non- 
linearity (a = — 1), it was impossible to identify two dif- 
ferent amplitude attractors and, thus, we were not able 
to construct dissipative DBs in this case. 

It is interesting to calculate the magnetization in a 
dissipative SRR array. In the direction perpendicular to 
the SRR planes the general relation 

B = (M) (H + M) (15) 
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FIG. 7: (Color online) Time evolution of £i n (T) (red-solid 
curve), compared with cos(57r) (black-dashed curve), and 
their sum (green-dotted curve), during two breather periods, 
for (a) a ID SRR array in the planar configuration at n = 15 
(A = —0.02, I — 3); (b) a ID SRR array in the planar config- 
uration at n — rib (A = —0.02, i = 3); (c) a ID SRR array in 
the axial configuration at n = 15 (A = 0.017, i = 2.4); (d) a 
ID SRR array in the axial configuration at n — rib (A = 0.017, 
t — 2.4). The other parameters as in Fig. 6. 

holds, with M. = SI /d 3 the magnetization (magnetic 
moment per unit cell volume) of the array. The earlier 
equation, with using Eq. ([3]), can be written in normal- 
ized form as 

B/B Q = cos(fir) +£i(T), (16) 

where B = e U c /Sui and I = ^ S 2 fl/e d 3 L. Eq. OH) 
may be used locally at each cell, with the three terms 
B/B , cos(fir), and £i(r) representing the instantaneous 
magnetic induction, applied magnetic field, and mag- 
netic response (local magnetization), respectively, in a 
specific cell. Negative magnetic response appears, as it 
is apparent from Eq. (|16p . whenever the second term on 
the right-hand-side of Eq. (|16[) larger in magnitude than 
the first one, and it has the opposite sign. In Fig. 7 
the time evolution is shown separately for each of the 
three terms of Eq. lfll)]) , i.e., the quantity ii n (r) (red- 
solid curve), cos(Ot) (black-dashed curve), and their sum 
B/B (green-dotted curve), for two different sites in both 
the planar and the axial geometries. Specifically, Figs. 
7(a) and 7(b) are for the planar geometry, for n — 15 (a 
site on the background, quite far away from the central 
DB site) and n = rib (DB central site), respectively, while 
Figs. 7(c) and 7(d) are for the axial geometry, for n = 15 
and n = rib, respectively. We observe that, in both ge- 
ometries, the SRR with low amplitude current oscillation 
(reduced nonlinearity) shows positive (paramagnetic) re- 
sponse (Figs. 7(a) and 7(c) for the planar and the axial 
geometry, respectively), while the SRR with high am- 
plitude current oscillation (enhanced nonlinearity) shows 
extreme diamagnetic (negative) response (Figs. 7(b) and 
7(d) for the planar and the axial geometry, respectively). 
Thus, in the breather or multibreather location the lat- 



tice has a negative magnetic response even though it is 
driven below resonance. That result can be extended 
to uniform solutions, where /„ = / for all n. Without 
nonlinearity such a solution always provides positive re- 
sponse below the resonance frequency (~ uji). However, 
the nonlinearity makes it possible to have two different 
coexisting and stable states (bistability) with low and 
high amplitude current oscillation, which are in phase 
and in anti-phase, respectively, with the applied mag- 
netic field. Thus, it is possible to obtain from an SRR 
array uniform negative response below resonance, by ex- 
citing all SRRs in the array to the high current amplitude 
state. These remarks are also valid for 2D SRR arrays 
discussed in the next section. 



IV. DISCRETE BREATHERS IN 
TWO-DIMENSIONAL SRR ARRAYS 

Although most of the methodology and techniques for 
DB construction has been developed for the ID case, 
there have been some studies of higher-dimensional sys- 
tems. Importantly, a rigorous proof of the existence of 
DBs in higher-dimensional nonlinear lattices was given 
in Ref. Q. Numerical studies of DBs have been pub- 
lished for several simple nonlinear lattices as for example 
2D Fermi-Pasta-Ulam chains [M S3, El SI, Josephson 
Junction arrays [50j . Klein-Gordon chains [51(, discrete 
nonlinear Schrodinger systems [Hj], and also for a 
Morse lattice [s3 |. Since most of the present MMs are 
fabricated in 2D technology, it is necessary to extend the 
study of magneto-inductive DBs in two dimensions. We 
see below that both the Hamiltonian and the dissipative 
DBs are not destroyed by increasing the dimensionality. 

The frequency spectrum of linear MI waves in the 2D 
system can be obtained as in the ID case, by substituting 
a plane wave of the form q n , m — A cos(K x n + K y m — fir) 
into the linearized equations Eqs. © and (TTUj) in the 
absence of losses and applied field (eo = 0, 7 = 0) 

tt K = [1 + 2 X x cos(k x ) + 2X y cos(Ky)] , (17) 

where k x and the (normalized) components of 

the wavevector in the x— and y— direction, respectively 
(— 7r < Ki < 7r, i — x, y). In this case, the group velocity 
is not, in general, in a direction opposite to the phase 
velocity [34]]. Notice that the bandwidth Af2 depends 
both on the magnitude of the coupling parameters X x 
and A y and their sign. Consider, for example, the case 
— A« = X (isotropic lattice in the planar geometry). 
Then Afl ~ 4|A| for |A| <C 1, which is larger than that of 
the corresponding ID system by a factor of two. Thus, in 
this case, the ID Hamiltonian DBs with frequencies very 
close to the ID band may not survive in the 2D case. 
Typical dispersion curves (i.e., contours of the frequency 
as a function of k x and K y ) of linear MI waves for isotropic 
2D SRR arrays in the planar geometry, anisotropic 2D 
SRR arrays in the planar geometry, and anisotropic 2D 
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SRR arrays in the planar-axial geometry, are shown in 
the left, middle, and right panels of Fig. 8, respectively. 




FIG. 8: Dispersion curves of the frequency spectra of lin- 
ear magnetoinductive waves (along with their corresponding 
density plots) on the k x — n y plane, for \ x — X y = —0.01 
(left panel); \ x = — 0.01, X y = —0.013 (middle panel); and 
A^ = —0.013, A H = 0.01 (right panel). In the density plots, 
darker color indicates higher fi K . 

For the construction of DBs in 2D arrays, we used 
the same methods as for the corresponding ID arrays. 
We again consider two different geometries of the SRR 
arrays, the planar geometry (Fig. 2(a)) and the planar- 
axial geometry (Fig. 2(b)). In the former geometry, the 
coupling parameters X x ,X y are both negative. However, 
they may differ in magnitude, i.e., \X X \ ^ \X y \ as a re- 
sult of unequal lattice constants d x and d y or resulting 
from the different orientation of the SRRs. In the latter 
geometry, the coupling parameters have opposite signs, 
i.e., X x < and X y > 0, which can be regarded as a case 
of generalized anisotropy. Their magnitude may however 
be different, depending again on the lattice constants d x 
and d y . In the following we are mainly concerned with 
single-site bright DBs for both geometries, nonlinearities 
(self- focusing and self-defocusing) , and several DB fre- 
quencies and pairs of coupling parameters X x , X y . How- 
ever, one may obtain many different types of DBs by just 
choosing the appropriate initial condition (trivial DB). 
The 2D array size used in the calculations is typically 
N x N = 15 x 15. 




FIG. 9: (Color online) Snapshots of two-dimensional Hamilto- 
nian single-site bright discrete breathers (taken at maximum 
amplitude) , for isotropic split-ring resonator arrays in the pla- 
nar geometry (X x — X y — A) with (a) a = +1, £l b = 0.952, 
A = -0.020; (b) a = +1, Q b = 0.938, A = -0.024; (c) 
a = +1, fl b = 0.881, A = -0.040; (d) a = -1, Q b = 1.082, 
A = -0.025; (e) a = -1, fl b = 1.036, A = -0.012; and (f) 
a = -1, Q b = 1.011, A = -0.004. 

Typical Hamiltonian DB profiles in 2D isotropic (X x — 
X y ) SRR arrays in the planar geometry are shown in 
Fig. 9 for both self-focusing (Figs. 9(a) - 9(c)) and self- 



defocusing (Figs. 9(d) - 9(f)) nonlinearities. The DB fre- 
quencies and coupling parameters are given, for each DB, 
in the caption of Fig. 9. Notice that for self-focusing non- 
linearity (a = +1), the DBs are staggered in both the x— 
and the y— directions, while for self-defocusing nonlincar- 
ity (a = —1), they are unstaggered in both the x— and 
the y— directions. We have constructed exact DBs also 
in the case of an array in the planar geometry with mod- 
erate anisotropy, i.e., Aa; ^ X yi for both nonlinearities 
(a = ±1). Typical DB profiles with the coupling param- 
eters differing by approximately 10% are shown in Figs. 
10(a) and 10(b) (upper panels) for a = +1 and a = — 1, 
respectively. The anisotropy does not change the stag- 
gered/unstaggered character of these DBs, which remain 
staggered (unstaggered) in both directions for a = +1 
(a = —1). The lower panels of Figs. 10(a) and 10(b) 
show the density plots of the corresponding DB profiles 
shown in the upper panels. In those density plots, where 
a darker color indicates a higher current amplitude, a 
DB which is staggered in both directions appears as a 
checkerboard rotated by 45 degrees with respect to the 
array. 
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FIG. 10: (Color online) Snapshots of two-dimensional Hamil- 
tonian single-site bright discrete breathers (taken at maxi- 
mum amplitude), for anisotropic split-ring resonator arrays 
in the planar geometry with (a) a = +1, Q b = 0.938, 
X x = -0.024; X y = -0.027; and (b) a = -1, fi 6 = 1.082, 
X x = -0.025, X y = -0.028. In both (a) and (b), the lower 
panel show the density plot of the discrete breather profile 
presented in the corresponding upper panel. In the density 
plots, darker color indicates higher breather amplitude. 

Typical Hamiltonian DB profiles in 2D anisotropic 
(I A x | 7^ I Ay |) SRR arrays in the planar- axial geometry 
are shown in Figs. 11(a) and 11(b) (upper panels) for 
a = +1 and for a = — 1, respectively, along with their 
corresponding density plots (lower panels). The cou- 
pling parameters differ in magnitude by approximately 
10%. In this case we observe that for a = +1 (a = — 1) 
the DB is staggered (unstaggered) along the x— direction 
(y— direction), while it is unstaggered (staggered) along 
the y— direction (x— direction). Thus, the change in ei- 
ther the sign of the nonlinearity a or the sign of the 



coupling parameter Xy, leads to a change in the stag- 
gered/unstaggered character of a DB in the y— direction. 
Specifically, by changing a from +1 to -1 (for X y > 0), 
a DB unstaggered in the y— direction becomes staggered 
in that direction, while by changing a from -1 to +1 
(for X y > 0), a DB staggered in the y— direction be- 
comes unstaggered in that direction. Also, by changing 
A j, from positive to negative (for a = +1), a DB un- 
staggered in the y— direction becomes staggered in that 
direction, while by changing X y from negative to positive 
(for a = +1), a DB staggered in the y— direction be- 
comes unstaggered in that direction. The linear stability 
of all the Hamiltonian DBs presented in this section is 
checked through Floquet analysis (finding the eigenval- 
ues of the Floquet matrix), and they were found to be 
linearly stable. 




FIG. 11: (Color online) Snapshots of two-dimensional Hamil- 
tonian single-site bright discrete breathers (taken at maxi- 
mum amplitude), for anisotropic split-ring resonator arrays 
in the planar-axial geometry with (a) a — +1, fij, = 0.938, 
\ x = -0.024; X v = 0.021;' and (b) a = -1, Q b = 1.082, 
X x = -0.025, X y = 0.022. In both (a) and (b), the lower 
panel show the density plot of the discrete breather profile 
presented in the corresponding upper panel. In the density 
plots, darker color indicates higher breather amplitude. 

Typical examples of dissipative DB profiles in 2D 
isotropic (A^ = X y ) SRR arrays in the planar geome- 
try, 2D anisotropic (X x ^ X y ) SRR arrays in the planar 
geometry, and 2D anisotropic (\X X \ ^ \X y \) SRR arrays 
in the planar-axial geometry, are shown in Figs. 12(a), 
12(b), and 12(c), respectively, for a = +1. These pro- 
files are actually snapshots at some specific instant dur- 
ing the DB motion. Here, just like in the ID case, both 
the background and the DB (i.e., the central DB site) 
oscillate with low and high current amplitudes, respec- 
tively, at the same frequency f2b = 17. The stability of 
the dissipative DBs has been checked by adding small 
perturbations of the order of 10 -2 and following their 
time evolution for long time intervals (over 10 3 Tf, time 
units). In all cases it was found that the DBs are not 
destroyed by the perturbation but, instead, they restore 
their unperturbed shape. 
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FIG. 12: (Color online) Snapshots of two-dimensional dissi- 
pative (single-site, bright) discrete breathers (taken at maxi- 
mum amplitude), for a = +1, Qt = 0.92, eo = 0.04, 7 = 0.01, 
6£ = 2 and (a) X x = X y — —0.02 (isotropic split-ring resonator 
arrays in the planar geometry); (b) X x — —0.02, X y = —0.023 
(anisotropic split-ring resonator arrays in the planar geom- 
etry); (c) X x = -0.0124, Ay = 0.0094 (anisotropic split-ring 
resonator arrays in the planar-axial geometry). 



The magnetic response of the dissipative 2D arrays can 
be calculated as in the ID case, with the help of Eq. 
(fT6)) applied locally at each array cell (n, to). In Fig. 
13 we plot separately the time evolution of each of the 
three terms of Eq. (fT6|) . i.e., the instantaneous magnetic 
induction, the applied magnetic field, and the magnetic 
response, at the central DB site (Figs. 13(b), 13(d) and 
13(f), high amplitude current oscillation) and the site 
with n, to = 3,5 (Figs. 13(a), 13(c) and 13(e), low am- 
plitude current oscillation), for the three DBs shown in 
Fig. 12. The results look the same for the first two 
cases, corresponding to isotropic 2D SRR arrays in the 
planar geometry (Figs. 13(a) and 13(b)), and anisotropic 
2D SRR arrays in the planar geometry (Figs. 13(c) and 
13(d)). They are also similar with those calculated for the 
corresponding ID case. The SRRs with low amplitude 
current oscillation ((n, to) = (3,5)) show positive (para- 
magnetic) response while the SRR with high amplitude 
current oscillation (central DB site) shows extreme dia- 
magnetic (negative) response. In the third case, however, 
for anisotropic 2D SRR arrays in the planar-axial geom- 
etry (Figs. 13(e) and 13(f)), the response of the SRR 
with high amplitude current oscillation is still diamag- 
netic but not negative. This is due to the very weak cou- 
pling which presumes relatively large separation of the 
SRRs, and thus a very low density SRR array. Guided 
from the previous results, we consider the possibility of 
constructing a region in a 2D SRR array with extreme 
diamagnetic (negative) response, surrounded by a para- 
magnetic background. For this purpose we may exploit 
a bright 2D multibreather consisting of a number of ad- 
jacent sites (SRRs) in the center of the array. For illus- 
tration, such a multibreather occupying a region of nine 
sites (including the central one) in an isotropic SRR ar- 
ray in the planar geometry, is shown in Fig. 14. We have 
checked that, indeed, the sites with high amplitude cur- 
rent oscillation show negative magnetic response, while 
the rest of them show positive magnetic response. Thus, 
it seems possible to create SRR-based MMs with distinct 
regions of opposite sign magnetic responses by exploiting 
multibreathers . 

The numerical results for both ID and 2D SRR arrays 
reveal that, at least for weak coupling, the amplitude and 
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FIG. 13: (Color online) Time evolution of £in(r) (red-solid 
curve), compared with cos(f2r) (black-dashed curve), and 
their sum (green-dotted curve), during two breather peri- 
ods, for a 2D SRR array (a) in the planar configuration at 
(n,m) = (3,5) (A* = X v = -0.02, I = 3); (b) in the planar 
configuration at the central breather site {Xx — X y = —0.02, 
I = 3); (c) in the planar configuration at (n,m) = (3,5) 
{Xx = -0.02, A y = -0.023, t = 2.6); (d) in the pla- 
nar configuration at the central breather site {X x = —0.02, 
X y = —0.023, I = 2.6); (e) in the planar-axial configuration 
at (n,m) = (3,5) (A* = -0.0124, \ y = 0.0094, I = 0.43); (f) 
in the planar-axial configuration at the central breather site 
{Xx = -0.0124, A H = 0.0094, £ = 0.43). The other parameters 
as in Fig. 12. 
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FIG. 14: (Color online) A snapshot of a two-dimensional dis- 
sipative multibreather (taken at maximum amplitude), con- 
structed for an split-ring resonator array in the planar geom- 
etry, for \ x = \ v = -0.02, uj b = 0.92, e = 0.04, 7 = 0.01, 
e e = 2, a = +1. 



the time-dependence of the central DB site and the back- 
ground are essentially those of the single damped-driven 
SRR oscillator. Only the first one or two sites neighbor- 
ing to the central DB site exhibit significant differences 
due to the coupling. The bistability of the single SRR 
oscillator, which is important for the construction of dis- 
sipative DBs and also for the creation of uniform states 
with either positive or negative magnetic response below 
resonance, is not restricted to frequencies close to reso- 
nance. In Fig. 15 we show the time evolution of i n of 



FIG. 15: (Color online) Time evolution of the normalized 
driving term cos(S7r) (black-dashed curve), and the current 
i(r) for the low and high current amplitude states (red-solid 
and green-dotted curves, respectively) for a single damped- 
driven SRR oscillator with 7 = 0.01, eg. = 2, a = +1 and (a) 

n = 0.8, £ = 0.2; (b) n = 0.5, £ = 1.2. 



two coexisting and stable states far from resonance, along 
with the normalized driving magnetic field. Dissipative 
DBs can be constructed by combining these two states 
into a trivial DB and continuing that solution for finite 
(non-zero) coupling parameters. Although for the cases 
shown the current of the high amplitude current state is 
rather large, where the saturation of the nonlinear term 
could be expected, Fig. 15 indicates that it is possible 
to obtain negative response below resonance by exciting 
the SRRs in that state. 



V. CONCLUSIONS 

We considered a simple model for nonlinear SRR-based 
MMs in one and two dimensions, where the nonlinear- 
ity arises from a Kerr-type dielectric which fills the SRR 
slits. Each SRR, which is modeled as a nonlinear RLC 
electrical circuit driven by an alternating voltage source, 
is weakly coupled to its nearest neighbors due to mag- 
netic dipole-dipole interactions through their mutual in- 
ductance M (magnetoinductive coupling). The sign of 
the coupling between neighboring SRRs depends on their 
relative orientation within the SRR array. 

We have constructed, using standard numerical meth- 
ods, many different types of Hamiltonian and dissipative 
DBs both in ID and 2D arrays for different nonlinearities 
(i.e., self- focusing and self-defocusing) , and different ge- 
ometries (planar and axial in ID, planar and planar-axial 
in 2D). We have also constructed DBs in 2D arrays with 
moderate anisotropy. Most of the DBs presented here 
are linearly stable under small perturbations. Dissipa- 
tive SRR arrays, driven by an applied magnetic field, of- 
fer the possibility to study their magnetic response with 
respect to that field. The induced current oscillations 
are proportional to the magnetic moments of the SRRs 
and, thus, to the local magnetization (magnetic response) 
of the array. We found that low (high) amplitude cur- 
rent oscillations are in phase (almost in anti-phase) with 
the applied field. Thus, depending on the array and the 
external field parameters, the magnetic response at the 
SRRs with high amplitude current oscillation can be neg- 
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ative. In this way, DBs can change locally the magnetic 
response of the array from paramagnetic to extreme dia- 
magnetic. By exploiting multibrcathers, it seems possi- 
ble to create SRR-based MMs with distinct regions of 
opposite sign magnetic responses. 

The bistability due to the nonlinearity is not restricted 
to frequencies close to resonance, but it is found to per- 
sist down to much lower frequencies. As a result, dissipa- 
tive DBs with frequencies and external field amplitudes 
in rather wide intervals can be constructed. Moreover, 
even far from resonance, the coexisting and stable low 
and high current amplitude states are in phase and in 
anti-phase, respectively, with the applied magnetic field. 



Thus, it seems possible to get uniform solutions at these 
low frequencies, which provide either positive or nega- 
tive magnetic response below resonance. This could be 
achieved by exciting all the SRRs in the low or high am- 
plitude state, respectively. 
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